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Multivariate analysis is that branch of statistics which is devoted to 
summarizing the relationships of sets of dependent variable;^. It Includes, 
for example, principal cdmpoient analysis or discriminant analysis, where 
the problem is developed in the statistical context of determining a linear 
combination of a given set of variables that has a larger variance than any 
other linear combination, or that best differentiates among several groups. 
The detailed theory and trends of extensive use particularly in behavioral 
science research will be introduced by such books as those by Cooley and 
Lohnes.[1971], Rulon [1967] and Tatsuoka [1971]. 

Those techniques will be useful in the realm of EEG pattern analysis, 
since the individual patterns may be assumed to be single points located in 
a multidimensional space. Hence, several multivariate techniques have al- 
ready been taken up recently in this field: step-wise discriminant anal'-sis 
[SWDA] to average evoked potentials [Donchin et al., 1970], component analy- 
sis to spectra, and the like. Particularly, it is vith the hope that appli- 
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cation of multJ.variate analysis will facilitate an dealing with taxonomic 
or classification problem of EEC. 

The purpose of this report is to describe sections of our studies in 
making use of the multivariate mathematical models to quantitative discrim- 
ination of EEG patterns among the groups of normal and several types of 
mentall retarded. We outline the preliminary approach in applying the 
techniques of principal component analysis and discriminant analysis. 

1. MATERIAL AND DATA PROCESSING 

Thirty-eight EEG pattern samples of 10 second epoch during resting 
coiidition were extracted from three groups of mentally retarded, namely, 
the predisposed [Group P] , the exogenous [Group E] and the Down's syndrome 
[Group D], who ranged in age frora 7 to 41 years. Those typical samples of 
each group, 11, 13 and 14 respectively, were selected from more than four 
hundred cases of the mentally retarded that were clinically examined. The 
subjects with evidence of any ei^ileptic and other neurological signs were 
excluded in selection. 

As controls, 32 norma?, samples ranging from 4 days to 20 years of age 
were also used [Group N] , which consisted of the younger 17 samples ranged 
up to 6 years [Gxoup Ni] and the other 15 samples matched in age to the re- 
tarded [Group N2]. These were extracted from the Gibbs' Atlas [1951]. 

Data processing was performed in three steps. In the beginning, EEG 
patterns were digitized by A/D conversion of sample waves on the magnetic 
tape and were punched automatically in 8 bit paper tape for computer pro- 
cessing, using TEAC R-400/ ATAC 501-10/ TH-SOO system at the laboratory. 

The binary coded decimal outputs thus obtained were then used to 



computer processing on multiple variables, listed in Table 1, to gain over- 
all infoimations of each pattern. Figure 1 shows the flov of data process- 
ing except that in auto/cross correlation and spectrum analysis, and its 
computer program is also given in the Appendix. 

The final step is to apply the component analysis and the discriminant 
analysis to sample values of variables obtained. Since the variables to 
use should be limited in number at these procedures, due to the limit of 
computer mq.mory, 26 variables were selected. Those are marked by Xj to X26« 

2 . COMPONENT ANALYSIS 

Suppose the random p-dimensional vector = [xi ,X2, . . ,Xp] has the 
variance and covariance matrix E. We shall assume that the mean vector is 
0 and x's have the unit variance. The object of component analysis is to 
economize in the number of variates, and for that, is to seek for a linear 
combination of type Z = a'A" which maximizes variance. 

Let a be a p-dimensional column vector such that a 'a = 1. Then the 
variance of Z is 

E[a'X]^ = E[a'XX'a] = a'la = a'Ra [1] 

where R is the correlation matrix. 

To determine the normalized linear combination a'X with maximum vari- 
ance, we must find a vector a satisfying a'a - I which maximizes [1]. In 
order to get a solution, we should seek for a satisfying 

[R - XI]a n 0 [2] 
where X is a Lagrange multiplier. If a satisfies [2] and a'a = 1, then the 



variance of a'X is X, Thus for the maximum variance we should use in [2] 
the largest x, namely, Xp 

Let ui be a normalized solution of 

[R - XjIJa « 0 . [3] 

•Then Zi « ai'X is a linear combination with maximum variance, and is called 
the first component. Furthermore, we may find another vector az correspond- 
ing to the second largest root X2 of [2], such that Z2 - az'X has maximum 
variance of all linear combinations uncorrected with Zj . Z2 is called the 
second component. This procedure is carried on, and we may thus transform 
to new variates Zi, Z2,.., which are uncorrected and have variances Xj, 
X2»..» Xp in decreasing order. ^ 

The results are shown in Table 2. It describes the coefficients of 
first seven linear combinations, which we obtained when applying the compo- 
nent analysis to 70 EEG samples simul" jously. The leading four extracted 
components, Zj, Z2, Z3 and Zi+ account respectively for ?3.0, 16.6, 13.3 and 
8.9 per cent of the total variance, and evidently 75 per cent are accounted 
for by the seven components given in Table 2. 

Multiplying each coefficient by /\j we have the correlation coefficient 
r^j- of ith variable and ^'th component; therefore, the signs of the coeffi- 
cients and their relative magnitudes are useful to examine the nature of com- 
ponents. 2 In the present results, it may be observed that Zi sums up infor- 
mation on "general development of EEG". Likewise, Z2, Z3, Z^,.. may be named 

^ To get a solution of [2] with a'a » 1 we. must have H - \I singular; in 
other words, X must satisfy |i? - xi| * 0. The function \R - XJ| is a poly- 
nomial in A of degree p. Therefore, the equation |i? - XI| = 0 has p roots; 
let these be Xj >_ X2 ^ . . . >^ Ap. 
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respectively "frequency versatility", "amount of irregular slowing", "a- 
mount of beta-activities at the occipital region", and so on [See Figure 2]. 

Figure 2 depicts the scatter diagrams representing relations between 
component scores and ages for Groups P. E, D and N. The lines and broken 
lines are showing the regression of the scatter diagrams for Group N. What 
is evident from these diagrams is that Groups E and D are retarded in regard 
to Z2 as compared with Group N. Furthermore, it is clear that Group D keeps 
high amount of occipital beta-waves, contrary to the result of decreasing of 
; those waves after 15 years in the other groups. However, definite tendencies 
cannot be observed in Group P; that will be one of topics for further dis- 
cussion. 

Thus, owing to reduction in the dimensions, the classification of EEG 
patterns may be discussed economically in terms of a set of fewer new var- 
iables, namely, components. But we think it is satisfactory to consider that 
the discriminant analysis will be more effective vo the classification prob- 
lem, which will be described in the next section. 

3. DISCRIMINANT ANALYSIS 

Suppose we have the the vector of p-dimensional measurements X' = [xj , 
X2i..» Xp] on an individual. : We shall now consider the assignment of that 



The covariance of and A is mciyx)X'] = a^^ViX'] = a^*R, The vari- 
ance cf is then the correlation coefficient = [ri^',... r^^] 
should be 

since Ra^ n \^a^ is derived from \R - \^l\a^^ 0, 



individual into one of two normal populations, namely, ttj : N[pi. z] and 
^2 : N[P2. I], where = m^] is the vector of means of the tth 

population [i = 1. 2] and I is the matrix of variances and covariances of 
each population. 

. In this case, if the observation X is actually from tt^. the linear 
combination Z - a'X should be distributed according to one-dimensional nor- 
mal distribution N[a\. a'la]. Then the problem is to classify into either 
^1 or .2 to minimize the distance \a'X - a'p^l , between \a'X - a'pij and 
\a'X - c\2\' 

Since the probability of misclassification, in this case, is the 
r.onic decreasing function of the Mahalanobis' distance between tii and t, 
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[a'ui - a'M2]^ 
a2 = 

a'Za ... 

the most appropriate weight a' = [a: a^] may be obtained by seeking for 

a so as to maximize under the restriction of 

a'lc = 1 . [5 J 

Thus Wr find 

a = KE-l[ui - [6] 
where K is the constant. The linear function 2 = a'X, thus obtained, is the 
well-known discriminant function, that is to differentiate best the obser- 
vations from two populations, iri and itj. 

The same result may be reached by a different route. The ith normal 
density function is 

^ exp[- l[X - ui] 'Z'HX - H)] [7] 

[2Tr]2 |j:|2 
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Then the ratio of densities is 

PlW exp[- i[X - yi]'E-l[^ - yi]] 

» f [8] 

P2W exp[- ^[X - \i2]'^-Hx - ViU 

« exp[;f'E-i[yi - \i2] - j[yi ♦ y2] £-Uyi - ^2]]- 

The region of classification into wj, which we denote by Ri, is the set of 
for which [8] is i k [for k suitably chosen]. Since the logarithmic function 
is monotonic increasing, the inequality can be written in terms of the log- 
arithm of [8] as 

X'Z'Hvi ' Ml] - j[yi + yz] Z'Hvi - y2] > log k . [9] 

If we denote the left-hand side as U, and if has the density [7] 
[i - I, 2], che best regions of classification are given by 

Ri: U ^ log k 
Ra: U < log k . 

If a priori probabilities qi and q2 are known, then k is given by 

q2 C[l/2] 

k« ^ [11] 

qi C[2/l] , 

where C[l/2] is the cost of misclassifying an individual from ^2 as from 
ffi, and C[2/i] is that in the opposite direction. 

In the particular case of the two populations being equally likely 
and the costs being equal, k = 1 and log k = 0. Then the regions of classi- 
fication into ffj and ^2 are respectively 

Ri: U 2 0 

R2J U < 0 . f^^J 
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If the first term of U will be denoted as a. namely, 

a = Z"^[Mi - M2] . [13] 

then we obtain 

U = ^'a . + ii2J 'a . [14] 

The first term is the discriminant function previously introduced. 
Transforming U sligiitly, we get 

u= [^'E-ipi -^p'lZ-ipi] - [^'z-ipa - i.p2rip2]. [15] 

Let us now define 

H ^ X'Z'\ - jV'iZ-hi . [16] 

then U may be stated as 

U = fi - f2 . [17]3 
We may then rewrite the classification proceaure as 



Ri: fl > f2 

fl < f2 . f^^J 



With using f^'s. we may proceed to the classification among m popula- 
tions [m > 3]. The population Ttj^ corresponding to the greatest f^/. say f^, 
among m f^'s will be the most appropriate population, to which X should be 
allotted. In sich a case, the probability of correct classification into 
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exp[f^ - max f^] 
Pi = 

Z exp[f^ - max fA . 



3 The Mahalanobis'. distance between X and the centroid of 7t^, if we denote 
this by D^^, should be d^i = [x - p^] 'r 4^ - v>:J - X'Z'^X . 2f^. Hence, we 
may also write f^ « ^[X'Z-'^X - D^^] and U « fi . ^2 = ji^h - d2i]. 



Things to be investigated will be divided into two cases: the one is 
allotting an individual EEG sample to one of two populations, namely, the 
normal and the mentally retarded in general, and. -the another is allotting 
to one of four populations: P. E. D and N. These are schematically illus- 
trated in Figure 3. 

The purpose of the analysis is to seek for p linear functions of the 
variables, f^. i equals one to p. so that a -sample_Qhs|rvation can be al- 
lotted to appropriate one of p populations, according to which of the f 's 
is the greatest when the sample values are substituted. Therefore, f^ 
might be called a measure of proximity to population ti^. As is evident 
from the upper diagram, in case of two populations, where to allocate a 
sample will be decided according to the value of the function U which we 
have by taking from fp If the value of U is positive, the sample 
should be allotted to the mentally retarded, and if negative, to the nor- 



mal . 



The table in Figure 4 shows the weights and the constant term of U. 
which we found by computation in case of two populations, using the sample 
data. At tha . Group N2 alone were used for the normal, for the sake of 
matching in age to Groups P. E and D. The values gWen to the whole samples 
by the discriminant function Z obtained in this way are also distributed in 
Figure 4. It is evident that two distributions for the mentally retarded 
samples and fjr the normal are clearly separated. Therefore, the probabil- 
ity of misclassification seems to be estimated as exi:remely low. 

When we proceed to discrimination among four populations, however, re- 
sults are more complicated. Table 3 shows the coefficients and the constant 
terms of the four linear functions fi. f2, fg and f^. which we obtained by 
computation. From the values given to the samples by four functions, prob- 



^ ton AVAIUBLE 



abilities of assignment to each one of four populations may be computed 
for each sample* These probabilities are tabulated in Table 4, The samples 
havihg a some amount of probabilities^ large and small i to be assigned tc 
the other groups are added by dashed lines • Probabilities of assignment to 
the normal are 1 for all normal EEG samples » but it can be seen that two 



samples of Group P and Group D, namely, P-7 and D-14, are misclassified to 
the popupation E with the probabilities of 0.618 and 0.945 respectively. 
Besides, complicated problem on classification may be pointed out for the 
sample E-6. The EEG patterns of those complicated samples classified to E 
are shown in Figure 6, comparing with E-8 that is typical of Group E. 

Thus, it is concluded that members o/ any groups of the mentally re- 
tarded were not misclassified to Group N at least in this study« However, 



can be drawn merely from sampling bias, which should be the subject for a 
future study. We assume that the sampling bias of the normal EEG was the 
primary factor affecting the result. 

Figure 5 gives the two-dimensional chart for fj through f3 with respect 
to which the individuals of the mentally retarded can be classified into 
three regions such that 



where U12 = - ^2* ^13 = fj - fa and U23 = f2 - fa* The space is divided 
by three boundary lines, U12 - 0> U13 =0, U23 • U13 - U12 = Oi intersecting 
at a single point. 

The dots, squares, and circles represent members of Groups P, E and D, 



the general veracity of this conclusion! is 



doubtful because such a result 




Ui3 > 0; U12 > 0 
U12 < 0; U23 > 0 
Ui3 < 0; U23 < 0 , 



[20] 
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respectively, and those added by sample numbers are the practically or 
probably misclassified ones given in Figure 6. Iv-hat is evident from Figure 
5 is that Group E lies in :lose proximity to Group D, when compared with 
the relation! of P to E and P to D. It is consistent with theresults of 
the component analysis shown in Figure 2, and such a result can be expected 
on pathological and empirical grounds [Hirai and Izawa, 1964], 

4. DISCUSSIONS 

t 

The following points are left as future problems: the one is what sort 
of multidimensional, variables should be introduced to identify and differ- 
entiate an individual EEG pattern exactly, and the unbiased sampling also 
should carefully be considered; that is another point. 

The results of multivariate statistical analysis may be said to depend 
finally on those two points. It may be true that the variables we intro- 
duced are mere preliminary ones; for that reason, further strict discussion, 
from physiological as well as statistical point, will be required on select- 
ing appropriate variables. As to the sampling, as well, it becomes a seri- 
ous problem that we used the normal EEG samples of the Gibbs' Atlas; those 
samples seem to be biased to fewer amount of fast waves and versatility, 
which will act in favor of discrimination from samples of the mentally re- 
tarded . 

The discriminant analysis may also be accomplished by finding the weights 
a that maximize the discriminant criterion X, defined by the ratio of the be- 
tween-groups to within-groups sums -of- squares of a linear combination a'X. 
The criterion X should be ' 

SS^, a *Ba 

T. TT-^^ [21] 



where B and V are the between-groups t-nd within-groups SSCP matrices, re- 
spectively, The necessary .condition for maximizing X reduces to 

[B - \W]o = C [22] 

which is equivalent to 

[r^B - \I]a s 0 [23] 

provided, as will generally be true, that V is non-singular. 

Thus, when this equation is solved, we get non-zero eigenvalues, which 
will be denoted as Xj, X2»<<» descending order of magnitude, and r as 

sociated eigenvectors ai, az,.-, a^- The elements of those eigenvectors 
may be used as combining weights to form r uncorrelated discriminant func- 
tions, the entire set of which constitutes the discriminant space [Rulon, 
1967; Tatsuoka, 1971] . 

For this attempt, the severe restrictions of normality and identical 
dispersion matrix in each group are not required; in this respect, this 
method exceeds that we used, apart from the discussions of discriminatory 
power. It seems to be a worthwhile subject to seek relationship between 
two methods in applying to EHG patterns. 
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Table 1. 

Variables employed for data processing. 



A: Intra-area Variables [LP § LO] 

LP ' LO 

mean wave frequency 
SD of wave frequency 
mean wave amplitude 
SD of wave amplitude 
weighted mean of frequency 
weighted mean of amplitude 
sum of D*A 

mode of frequency [1] 
mode of frequency [2] 
sum of 6 amplitude 
sum of a amplitude 
sum of Si amplitude 
total amount of 6 waves 
total amount of a waves 
total amount of 6i waves 
auto-correlation 
auto- spectrum 
power of 4-7 Hz 
power of 8 Hz 
power of 9-10 Hz 
power of 11-12 Hz 
power of 13-19 Hz 
peak frequency 



• B! Inter-area Va riables [LP-LO] 

35 cross-correlation 36 cross- spectrum 



^ using amplitude as weight 

^ using wave duration as weight 

^ D*A: wave duration multiplied by amplitude 

[1] frequency which shows the highest peak of sum of amplitude 
[2] frequency which shows the highest peak of wave numbers 

5 band-width: 9[4-8 Hz], a[8-13 Hz] § 6i [13-20 Hz] 
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Figure 1. Data processing of way.As 
and its flow-diagram 
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Table 2. 

Results of component analysis applied to EEG samples. 
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Figure 3. Schematic diagrams of the discriminant analysis 
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0.10985 


25 


-0.07969 


0.20190 


-0.07484 


26 


-0.13377 


0.33889 


-0.10360 



constant term of U -0.18768 
Mahalanobis* distance 44.11321 

prob. of mis classification 0.00000 
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Figure 4. Results of discriminant analysis 
in assignment into two populations 
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Table 3. 

Coefficients of four linear functions of the variates f. 
when a sample will be allotted to 
one of four populations. 





f, 
1 


f 


'3 


f 
^4 


1 


-3.21223 


-1. 70721 


-4.73777 




2 


5.51442 


1.50793 


5.75684 


-10 72541 


3 


3.1579K 


3.56681 


4.56734 




4 


1.87781 


-5.30259 


-3.32020 


6 '51740 


5 


-1.88480 


-0.36555 


-4 13029 




6 


-1.16677 


-1.0983*i 


-0.27543 


2 nf>444 


7 


-1.34622 


-0.55439 


1.15727 


0 ^87^0 


8 


-1.08231 


-2.64588 


0.44676 




9 


0.48790 


-0.28345 


-0.69710 

V . V ^ r Aw 


0 S3852 


10 


•0.31535 


-2.60035 


-1 .62966 


4 00583 


11 


-1.346C8 


-1.40560 


-1.97901 


4 05245 


12 


0.69760 


3.79400 


1 .90643 


-5 57898 


13 


-0.32422 


2.90533 


0.74993 


-2 98014 


14 


1.45954 


-2.00025 


-1.09966 


1.68953 


lb 


1 . 32677 


8.80853 


-3.23178 


-5.59030 


16 


0.20775 


1.60024 


4.44757 


-5.68913 


17 


1.35211 


-1.14910 


-1.28686 


1.20499 


18 


2.31651 


-0.27492 


-2.60211 


.0.96760 


19 


1.86774 


2.47589 


-3.25302 


-0.47992 


20 


-0.15522 


0.44112 


1.66944 


-1.82678 


21 


-0.08244 


-3.60310 


0.27595 


2.92482 


22 


-1.46977 


0.65849 


-1.00890 


1.44881 


23 


-0.82224 


1.83958 


-0.50786 


-0.51727 


24 


3.50552 . 


1.92631 


0.07033 


-4.30573 


25 


-1.09734 


-2.54092 


0.04611 


2.96390 


26 


0.45734 


-4.19572 


-0.63146 


3.39018 




-5.22886 


-3.81538 


-2.53466 


-11.07762 



The second term of function f. = X'I'^m, - ^ y.'E"% 
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Table 4. 

Probabilities of classification when samples will be 
allotted to each one of four populations. 







^ 




"3 


•"^ 


p 


1 


0.99993 


0,00001 


0.00007 


0.00000 




2 


0.99939 


0.00042 


0.00019 


0.00000 




3 


1.00000 


0.00000 


0.00000 


0.00000 




4 


0 99824 




U. UUiSi 


0.00000 




s 


0 99918 




U. 00059 


0.00000 




6 


0.97467 


0 00087 




0. 00000 




7 -^E^ 


0, 00876 


0 617Q1 


n T 7 0 c n 


n nnno9 




8 


1.00000 


0 00000 


U . OOOUO 


O.OOOOJ 




9 


0.80300 

W t W WW w w 


0 Ol2f%d 

w . 1/ X 6 


n 1 QA 


n /\/\/\/\/\ 




10 


0,83053 

w g WWWWW 


0 OlOfid 


n 1 ^flflT 


n nnnnn 




11 


0.99989 


0 00007 


n nnni i 
U . 0001 1 


0. 00000 


E 


1 


0. 00001 

wa wwvA 




0. 000o4 


0.00000 




2 


0 00115 




0* 02o77 


0.00000 




3 


0 G0065 


0 qqc;qq 


0. OOooo 


0.00000 




4 


0.16471 


0.82509 


0.01020 


0.00000 




S 


0.00008 


0.99881 


0.00112 


0.00000 




6 


0.37708 


0.41097 


0.21195 


0.00000 




7 


0.00014 


0.79351 


0.20636 


0.00000 




8 


U. 00004 

W • WWWW*T 


0 0^009 


n n/^naA 
U. U0U}j4 


0. 00000 




9 


0.00007 


0.96553 


0.03440 


0.00000 




10 


0.00001 


0.97994 


0.02005 


0.00000 




11 


0.00006 


0.99872 


0.00121 


0.00000 




12 


0 . 00005 


0.93809 


0.06186 


0.00000 




13 


0 00090 

w • w ww7 w 


U . 7743D4 


0. 0U546 


0. 00000 


D 


1 


0.29930 


0.00114 


0.69956 






2 


0. 00526 


n 0009A 




0. 00000 




3 


0.00033 

W t WWWww 


0 onoK;A 


n OQQ1 1 
0 . 7771 1 


U. 00000 




4 


0. 14899 


0 0^7Q1 


n Q 1 T 1 n 


n nnnnn 




5 


0.00034 


0.02543 


0.97423 


0.00000 




6 


0.00010 


0.03347 


0.96643 


0.00000 




7 


0^00019 


0.18173 


0.81808 


0.00000 




8 


0.00202 


0.00227 


0.99571 


0.00000 




9 


0.00034 


0..?8947 


0.71019 


0.00000 




10 


0.01906 


0.00422 


0.97652 


0.00021 




11 


0.00000 


0.00026 


0.99974 


0.00000 




12 


0.00081 


0.00575 


0.99344 


0.00001 




13 


0.00023 


0.07799 


0.92179 


0.00000 




14 H-E 


0.00027 


0.94497 


0.0S476 


0.00000 



tnisclassification 
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Table 4. 

Probabilities of classification when samples will be 
allotted to each one of four populations. 

[continued] 







*^2 


P3 


P4 


N2 1 


0.00000 


0.00000 


0.00000 


1.00000 


2 


0.00000 


0.00000 


0.00000 


1.00000 


3 


0.00000 


0.00000 


0.00000 


1.00000 


4 


0.00000 


0.00000 : 


0.00000 


1.00000 


5 


0.00000 


0.00000 


0.00000 


1.00000 


6 


0.00000 


0.00000 


0.00000 


1.00000 


7 


0.00000 


0.00000 


0.00000 


1.00000 


8 


0.00000 


0.00000 


0.00000 


1.00000 


9 


0.00000 


0.00000 


0.00000 


1.00000 


10 


0.00000 


0.00000 


0.00000 


1.00000 


11 


0.00000 


0.00000 


0.00000 


1.00000 


12 


0.00000 


0.00000 


0.00000 


1.00000 


13 


0.00000 


0.00019 


0.00000 


0.999-81 


14 


0.00000 


0.00000 


0.00000 


1.00000 


15 


0.00000 


0.00000 


0.00000 


1.00000 




' ' ' ' » ' ' ' I ' ' ' ' ' ' ■ 

16 12 8 4 0 -4 -8 -12 



Figure S. The regions separating the three 
groups of the mentally retarded 
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C EH3 ANALYSIS ( HISTOGRAM METIHOD ) 

DIMENSION FNM(10,3) ,D(16) ,X(700,2) ,R(3) ,AV(2) ,SD(2) ,ASP(30,2) ,FRBQ 
1 (2) ,AMP(2) ,SP(2) ,SA(2) ,1H(2) ,BE(2) , 

2N(2) ,NUM(700) ,NSP(30,2) ,NAV(2) ,NSD(2) ,NIH(2) ,NAL(2) ,NBE(2) 
DO 10 J«l,3 

READ(2,101) (m(I,J),I»l,10) 

101 FORMAT (10A8) 
10 CONTINUE 

D=l 
J«l 

16 I1«0 

14 READ(2,102) (D(I) ,I»1,16) 

102 FORMAT (16F5.0) 
DO 12 I"l,16 

IF(D(I).EC.0.01) GO TO 13 
Il-Il+l 

12 X(I1,J)=D(I) 
GO TO 14 

13 I2«I+1 
N(J)«I1 

IF(J.BQ.2) GO TO 15 
J»2 

IF(I.EQ.16) GO TO 16 
I1«0 

DO 17 1-12,16 
Il-Il+l 

17 X(Il,J)iO(I) 

GO TO 14 1 

15 IF(N(1).I£.N(2)) NSMN(l) 
IF(N(l).Gr.N(2)) NS^(2) 
OOV-0.0 

DO 18 Jb1,2 

AV(J)«0.0 

SD(J)«0.0 

DO 18 >1,NS 

AV(J)«AV(J)+X(I,J) 

18 Sb(J)«SD(J)+X(I,J)**2 
DO 19 IxlfNS 

19 CXV<CV+X(I,1)*X(I,2) 
SN^S 

R (1) « (OOV*SN-AV (1) W (2 ) ) /SOTT ( (SD (1) *SN-AV (1) **2 ) * 
1(SD(2)*SN-AV(2)**2)) 
DO 20 J«l,2 
NOO 

FREQ(J)«0.0 

AMP(J)bO.O 

SF(J)«0.0 

SA(J)«0.0 

MrH(J)«0 

NAL(J)-0 

NBE(J)«0 

TO (J) -0.0 

AL(J)«0.0 
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BE(J)«=0.0 
DO 21 I«l,30 
ASP(I,J)bO.O 

21 NSP(I,J)*0 
DO 22 1=1,700 

22 NUM(I)bI 
32 I»l 

K«l 

25 I1«I+1 

IP(X(I,J) .LT,X(I1,J)) GO TO 23 
I-I+l 

IF(I.EQ.NS) GO TO 24 
GO 'SO 25 

23 X(K,J)-X(I,J) 
. XMIN1»X(I,J) 

KS»NUM(I) 

27 I«I+1 
IF(I.BQ.NS) GO TO 24 
11=1+1 

IP(X(I,J).GT.X(I1,J)) GO TO 26 
GO TO 27 

26 XMAX»X(I,J) 
K1«NUM(I) 

29 

IF(I.SQ.NS) GO TO 24 
11=1+1 

IF(X(I,J).LT.X(I1,J)) GO TO 28 
GO TO 29 

28 XMIN2=X(I,J) 
K2=NUM(I) 
Gl«FLQAT(rd-KS) 
G2"FLQAT(K2-KS) 

AfXMAX-XMINI-GI* (XMIN2-XMIN1) /G2 

Gl=60.0/G2 

MF=<;1 

FREQ(J)=FREQ(J)-K31 
AMP(J)=itflP(J)+A 
. SF(J)=SF(J)-K31**2 
SA(J)=SA(J)+A**2 
NO=NOfl 

IP(MF.17r.l) GO TO 30 
IF(MF.GT.30) GO TO 30 
ASP(MF,J)«ASP(MP,J)+A 
NSP(MF,J)«NSP(MF,J)+1 

30 NUM(K)=KS 
K-K+1 
JJJM(K)«K2 
XMIN1=XMIN2 
X(K,J)<"XMIN2 
KS=K2 

GO TO 27 

24 IF(K.LT.3) GO TO 31 
NS=K 

GO TO 32 

31 StMX) 
PREQ(J)=PREQ(J)/BN 
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AMP(J)«AMP(J)/BN 

SF (J) =SF (J) /BN-FREQ (J) **2 

SA (J) =SA (J) /BN-AMP ( J) '**2 

DO 33 1-4,7 

TH(J)«TH(J)+ASP(I,J) 

33 NTH(J)«Nni(J)+NSP(I,J) «{ 
DO 34 1-8,12 
AL(J)«AL(J)+ASP(I,J) 

34 NAL(J)»NAL(J)-»-NSP(I,J) 
DO 35 I«.13,19 
BE(J)-BE(J)-»>ASP(I,J) 

35 NBE(J)«NBE(J)-»>NSP(I,J) 
20 CONTINUE 

. DO 36 J«l,2 
AV(J)«0.0 
SD(J)«>0.0 
NAV(J)«0 
NSD(J)-0 
DO 36 I«l,30 
AV(J)»AV(J)+ASP(I,J) 
SD (J) -SD (J) +ASP ( I , J) **2 
NAV ( J) »ittV ( J) 4^SP (I , J) 

36 NSD(J)«NSD(J)-HgSP(I,J)**2 
COV-0.0 

NO»0 

DO 37 I'-l,30 

CCV-CCV+ASP (1,1) *ASP (1,2) 

37 N0bNCHNSP(I,1)*NSP(I,2) 

R (2) ■ (COV*30 . 0-AV (1) *AV (2) ) /SQRT ( (SD(1) *30 . 0-AV (1) **2) 
1*(SD(2)*30.0-AV(2)**2)) 

COV-NO 

DO 9 1-1,2 

AV(I)-NAV(I) 
9 SD(I)«NSD(I) 

R (3) - (CXJV*30 . 0-AV (1) *AV (2) ) /SQRT ( (SD (1) *30 . 0-AV (1) **2) 
1* (SD(2)*30.0-AV(2)**2)) 

WRr]:E(3,300) (FI»4 ( I, L), 1-1,10) 

300 FOFMAT(1H7,10A8/1H2,2X,IOHMEAM FREQ. »10X,2HSD,3X,9IMEAN AMP., 
110X,2HSD,2 (7X,5Hnim,7X,5HALPHA,8X,4HBEI!A) ) 

DO 38 J-1,2 

WRITE (3, 301) FREQ(J) ,Sr (J) ,AMP(J) ,SA(J) ,<ZH(J) , 
lAL(J) ,BE(J) ,NrH(J) ,NAL(J) ,NBE(J) 

301 F06MAT(1H2,7F12.4,3I12) 

38 CONTINUE 

WRITE (3, 302) (R (I), 1-1,3) 

302 F0PMAT(1H3,1X,11HC0RRELAT10N,3F12.4) 
IF(L.EQ.3) GO TO 11 

L-Lfl 
J-1 

IF(I2.EQ.17) GO TO 16 
. Il-O 

DO 39 1-12,16 
Il-Il+l 

39 X(r.,J)-D(I) 
GO TO 14 

11 STOP 
END 



- 26 - 



